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We propose an efficient numerical algorithm to solve Bogoliubov de Gennes equations 
self-consistently for inhomogeneous superconducting systems with a reformulated poly- 
nomial expansion scheme. This proposed method is applied to typical issues such as a 
vortex under randomly distributed impurities and a normal conducting junction sand- 
wiched between superconductors. With various technical remarks, we show that its 
efficiency becomes remarkable in large-scale parallel performance. 
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1. Introduction 

The mean-field approach is one of the most convenient and efficient ways to clarify 
relationship between microscopic descriptions and macroscopic phenomena in condensed 
matter physics. So far, this approach has been applied to various fermionic interacting 
many-body systems. Recently, inhomogeneous superconducting systems such as nano- 
scale superconducting wires and artificial/intrinsic Josephson junctions have been its 
intriguing application examples. 1-4 ^ The Bogoliubov-de Gennes (BdG) equation 5 -* is a 
theoretical starting-point in these systems. Its self-consistent solution gives information 
such as non-trivial quasi-particle excitation-spectra and inhomogeneous superconduct- 
ing gap. However, such a mean-field BdG approach has been regarded to be not practical 
for a long time, since its calculations require huge computation resources beyond the 
contemporary standard level of present computers. 

Very recently, a few groups 6-8 ^ have proposed a highly efficient numerical method to 
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solve the BdG equations by using the kernel polynomial expansion. 9 ) In these papers, 
the key idea is to expand Green's function with a set of the Chebyshev polynomials. 
The polynomial expansion drastically reduces the computation cost and has an excellent 
parallel efficiency in contrast to the conventional way, i.e., direct full diagonalization of 
the BdG-Hamiltonian. In the history of condensed matter physics, such a Chebyshev- 
based expansion of Green's function goes back to seminal studies by Tanaka et al. 10,n ) 
These authors calculated Green's functions for a generic Hamiltonian with use of various 
kinds of orthonormal polynomials. Sota et al. developed this idea as an oscillation- 
free Fourier expansion scheme. Afterwards, several papers have been published in terms 
of the polynomials expansion of Green's function. 9 ' 13 ) Thus, the polynomial-based ex- 
pansion has been an attractive numerical method in large-scale fermionic many-body 
systems. Specifically, the application to inhomogeneous superconducting systems is an 
important target since spatial profile of the superconducting gap has key information 
to predict their physical properties. 

The aim of this paper is to develop a fast and tractable method to self-consistently 
calculate the BdG equations for various inhomogeneous superconducting systems. We 
reformulate the polynomial expansion scheme in a more comprehensible manner and 
perform first self-consistent calculations on typical inhomogeneous superconducting sys- 
tems. At first, we claim that the polynomial expansion is applied to not Green's function 
itself but the spectral density of Green's function. By expanding the Dirac's delta func- 
tion on the spectral density, we obtain a full mean-field calculation scheme. The present 
tool, spectral- density polynomial expansion allows a straightforward numerical calcula- 
tion with the mean-fields, i.e., superconducting gap and/or general Hartree-Fock terms. 
Next, we claim that the mean- field itself converges faster than the local density of states. 
We focus on s-wave superconductivity to demonstrate its efficiency. We remark that our 
approach can be effective for more general cases including d-wave superconductivity in 
magnetic fields 14 ) and multi-orbital superconductivity with spin-orbit coupling. 

This paper is organized as follows. The polynomial-expansion reformula is given in 
Sec. II. The present style is theoretically more complete and easier to do programing 
on parallel cluster machines. In Sec. III., we demonstrate three typical examples; an 
s-wave disordered superconductor in the magnetic field, i.e., a vortex formation under 
randomly distributed impurities, a nano-scale superconductor-normal-superconductor 
(SNS) junction, i.e., superconducting proximity effects, and a numerical challenge, i.e., 
a reproduction of full temperature dependence of the gap amplitude from zero to T c . 
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Some convenient technical remarks for practical simulations are given in Sec. IV. Section 
V is devoted to the summary. 

2. Formulation 

We give a full formula based on the orthonormal-polynomial expansion to solve the 
BdG equations. The essence is that the spectral density of Green's functions can be 
expanded by orthonormal polynomials. Afterwards, we show the explicit expressions 
with use of the Chebyshev polynomials. The present formalism can be applied to any 
fermionic quadratic Hamiltonian with mean-fields. 



2. 1 Hamiltonian 

We start with the Hamiltonian associated with the BdG equations. Covaci et al. 6 ^ 
proposed a way to calculate the eigenvalues and the eigenvectors of this Hamiltonian 
without full diagonalization. We describe their formula in a more general way. 

Let us consider a Hamiltonian for a fermion system given as H = ^WS> /2. The 
column vector \1/ is composed of N fermionic annihilation q and creation operators 
c\ (i — 1,2,..., N), — ({ Cl },{c\}) T , where {cj = (c u c 2 , . . . , c N ) T and {cj} = 
(c\, 4, . . . , cjv) T . The row vector ^ f is also defined as ^ f = ({cj} T , {q} t ). The symbol T 
means transposition. The fermionic canonical anti-commutation relation leads [q, = 
5ij. The subscription i in q or cj indicates a quantum index depending on spatial site, 
spin, orbital, etc. The "Hamiltonian" matrix % is a 2iV x 2N Hermite matrix given as 

H { ' U \. (1) 
I B* -.4 T J 

where A and B are N x N complex matrices. These matrices have the relation given as 

it = a, B T = -B, (2) 

because of the hermitian property of H and the fermionic canonical anti-commutation 
relations. When we consider a superconductor, % corresponds to the mean-field 
Bardeen-Cooper-Schrieffer (BCS) Hamiltonian and B contains the superconducting 
gap. 
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2.2 BdG equations 

The BdG equations are regarded as th eigenvalue equation with respect to H ex- 
pressed as 

HfM = e y fw (7 = 1,2,..., 27V), (3a) 

(3b) 

The column vectors tt( 7 ) and vm are ^-component complex vectors. To solve the BdG 
equations is equivalent to diagonalization of % with a unitary matrix U (Ref. 15), 

U^nU = D, D = dmg(e 1 ,e 2 ,...,e 2N ). (4) 

The eigenvalues e c 's are not independent of each other, i.e., = — e i+N (i — 1, 2, . . . , N). 
The matrix elements of U lead as 

C/ i7 = U( 7 ) ; j t/i+Ar 7 = W( 7 ) 5 j. (5) 

£.3 Spectral density 

Here, we concentrate on a key quantity, i.e., spectral density (or discontinuity) d(ou), 
which is a 2N x 2N matrix, to solve the BdG equations using orthonormal polynomials. 
All essential physical observables are described by bilinear forms with respect to d(u). 

Now, let us define the Green's function as G(z) = (z — "H) _1 , which is a 2N x 2N 
complex matrix. With the use of the unitary matrix U, each component of G(z) is 
expressed as 

2N 

z — e-, 

7=1 ' 

If we set z = iu) n with the Matsubara frequency u n = (2n + l)//3, the above formula 
corresponds to the temperature Green's function. The retarded and the advanced 
Green's functions are, respectively, defined as 

G R {u) = lim G{u + ir]), (7a) 

G A {u) = lim G(u-ir]). (7b) 

The spectral density 16 ) is given as a difference between the retarded and the advanced 
Green's functions, d(u) = G R (u) — G a (uj), whose matrix elements are expressed as 

27V 

d(uj\ = -2vri V U ai mj(uj - e 7 ) . (8) 



2N 1 

G afi (z)=Y,U<ryU*p y —- (l<a,P<2N). (6) 
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In order to obtain physical observables (e.g., density of states) from d(u>), we introduce 
the following useful 2iV-component unit-vectors e(i) and h(i) (1 < % < N), which are, 
respectively, defined as 

[e(i)] 7 = <J i>7 , [/i(i)] 7 = <W, 7 - (9) 

Specifically, employing e(i) and ^i(i), the elements of the column vectors tt 7 and w 7 are 
rewritten as 

«f 7)>i = [C/ t /i(i)] 7 - (11) 
Furthermore, the local density of states with respect to the site % is given as 

N(u,i) = -^e(i) T d(uj)e(i), (12) 
with use of d{u) and e(i), since N(uj,i) is defined as 

2N 

JV ( w > i ) = 5^l U (7),i| 2<J ( £<;-e 7)' (I 3 ) 
7=1 

TV V 

= 1] l u Ci),i| 2<J ( w - e i) + J] !%)^! 2<5 ( w + e i)' (I 4 ) 

A typical self-consistent BdG calculation for a superconductor requires two types of 

mean-fields (c\cj) and (qCj). These mean-fields can be expressed as, 

1 /"°° 

(c < tcj) = ~2^i 7 _ <WMeO') T d(a;)e(0, (15a) 
(ciCj) = ~— / duf(u)e(j) T d(u)h(i), (15b) 
with f(x) = 1/ (e^ + 1). Here, (3 is the inverse temperature j3 = 1/T. 

2.4 Orthonormal polynomial expansion 

In this paper, we focus on an orthonormal polynomial (f> n (x) with interval [-1,1] 
(n — 0, 1, . . .). In principle, various orthonormal polynomials are applicable to solve the 
BdG equations (as shown in Ref. 11). A set of polynomials 4> n (x), which is assumed to 
be a real function with respect to x, fulfills the relations 

W{x) 



5(x - x') = V ^^(^ n (i'), (16a) 

Wj n ,m = J (f)„(x)(f) m (x)W(x)dx. (16b) 



J. Phys. Soc. Jpn. 



FULL PAPERS 



A recurrence formula is generally given as 17 ) 

<Pn+i(x) = (a n + b n x)(j) n (x) - c n n _i(x). (16c) 

In order to confine the eigenvalue range inside the interval, we rescale the energy scale 
of % by the following manner 

*=*Z*f, ^ = hzl, (17) 
a a 

where a = (£ max - E min )/2 and b = (E nmx + E miQ )/2 with E min < e 7 < E m3uX . It seems 
that these relations require a tough computing task to obtain E mSuX and E min . However, 
their rough estimations are practically enough. We employ a convenient criterion to 
determine E max and E m i n from physical information. See Sec. 4.2 for detail. 

Now, let us derive a formula using the polynomial expansion. At first, we define a 
matrix form by polynomial functions as 

2N 

4>n{£)} =Vc/ Q7 t/^(g, (is) 

where 0n(£ 7 ) is well-defined in the interval £ 7 G [—1, 1]. This implies that w-integrals 
in Eqs. (15a) and (15b) are also bound in the finite energy range. According to a 
similar manner to Eq. (17), these integral intervals become [—1, 1], i.e. u = ax + b with 
x E [—1, 1]. Substituting the right hand side of Eq. (16a) for the definition of d(tu), we 
have 

p?d(u) q = -—jt ^U^P T q n , (19) 

n=0 n 

for arbitrary 2iV-component real vectors p and q. A sequence of the vector q n {= 
4> n (iC)q) is recursively generated by 

q n+1 = (a n + b n K)q n - c n q n _i (n > 2), (20a) 

qi = M£)q, <?o = 0o(/C)q. (20b) 

The coefficients of Eq. (20a) are the same as the ones in Eq. (16c). Accordingly, the 
mean-fields (15a) and (15b) are expressed as 

fe> = E e w Te ^)-' ( 21a ) 

n=0 Wn 

(c lCj ) = ^e(j) T ^W — , (21b) 

n=0 Wn 
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where 

T n = J dxf(ax + b)W(x)(P n (x), (22) 

e n (i) = <j> n (K.)e(i), h n (i) = 4> n {fc)h{i). (23) 

Here, it should be noted that T n does not depend on the index i. Therefore, the calcula- 
tion of T n can be done before any self-consistent calculations. The essential mathemat- 
ical relations are Eqs. (16a)- (16c). Many useful formulae about orthogonal polynomials 
applicable to physical problems are found, for example, in Refs. 18 ' 19 ^ 

Hereafter, as shown by by Covaci et al.f 1 we use the Chebyshev polynomials, 17 ^ i.e., 

<f>n(%) — cos[n arccos(x)], (24a) 

W(x) = -=L=, w n = ^(l + 5 n0 ), xe[-l,l]. (24b) 

The coefficients in the recursive formula (16c) are a n = 0, b n = 2, and c n = 1. The 
vector form of the formula associated with Eq. (20a) is given as 

q n+ i = 2Kq n - <?„__! (ra>2), (24c) 

with q = q and q 1 = ftq. At the zero temperature (j3 — > oo), we find that 

7o = 7r — arccos(— b/a), (25a) 

Tn^o = sin i narccos (~ 6 / a )] (25b) 

Th 

Equation (24c) allows the evaluation of the right hand side of Eq. (19) without any direct 
diagonalization of "H. Specifically, we find that the calculations of the mean-fields (21a) 
and (21b) do not require any heavy computation in contrast to matrix diagonalization. 

Finally, we remark that the approach shown in this section is applicable to solve the 
other equations in condensed matter physics such as the Kohn-Sham equation in the 
density functional theory with real-space formalism. 

3. Results: Numerical Demonstrations 

We demonstrate three examples of self-consistent calculations in inhomogeneous su- 
perconducting systems by using the above spectral-density polynomial expansion. In 
this section, we focus on a sigle-band s-wave superconductor. Therefore, the quantum 
index i(J) introduced in the previous section indicates a single spatial site. Then, the 
spin index f (|) is separately written. The label (i, j) simply means a spatial site in 
real 2D space. Together with the BdG equations, the gap equation for s-wave supercon- 
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ductivity is given as Ajj = (ci,4.Cj,t) ■ m numerical calculations, its right hand side is 
expressed by using Eq. (21b) with truncation in the infinite summation. The maximum 
number of the summation is written as n c . For simplicity, we consider a simple square 
lattice tight-binding model only with nearest neighbor hopping. The hopping magni- 
tude is t. In all simulations, we adopt n c = 1000. We find that this number is enough 
from convergent tendency except for the calculations of the temperature dependence of 
A ir 

3.1 Vortex lattice Solution 

The first example is a vortex solution and related issues. The electromagnetic re- 
sponse of a superconductor under the magnetic field has been examined on the basis of 
the BdG equations . 20 ~ 22 ) In type II superconductors, the low-lying vortex-core excita- 
tions can be studied by taking account of the presence of vortex lattice. 

Here, we investigate a two-dimensional N x x N y site system with periodic boundary 
condition. This system has a vortex square lattice. Moreover, we introduce randomly 
distributed impurities at the zero temperature. The parameters are set as follows: V^j = 
—2.2t5ij, the chemical potential fj, = — 1.5t, the system size N x x N y = 64 x 64, and 
the impurity potential VJ mp = t. We successfully calculate the gap amplitude as shown 
in Fig. 1. We note that a calculation with 20 times iterations takes about 40 minutes 




10 20 30 40 50 60 70 

(a) (b) 

Fig. 1. (Color online) (a) Spatial modulation of the gap amplitude in a two-dimensional vortex 
lattice system with ten impurities (b) Cross section of the spatial modulation at y = 32. 

by a desktop computer with 8 CPU cores (Intel Xeon X5550 2.67 GHz x 2). Since the 
Hamiltonian is a sparse matrix, the calculation needs little computational memories. 
The separate calculations of A^- is performed on each CPU core. The communication in 
this case, which is a one-to-all communication, is needed only when updating Ay. We 




confirm that these advantage points are still true in more general cases such as <i-wave 
and multi-orbital superconductors. 23 ' 

3.2 SNS Junction System 

The next example is the proximity effects in an SNS Josephson junction, in which 
the gap symmetry of the superconducting electrodes is assumed to be s-wave. Joseph- 
son junctions and superconducting weak links have attracted great attention in vari- 
ous research fields such as superconducting device engineering 24 ^ including Josephson 
qubits, 25 ) superconducting transport studies in superconducting wires, 2&> fundamental 
studies on unconventional nano-superconductors, 7 ' 27-29 ^ and so on. Numerical simula- 
tions for the BdG equations has been often carried out to study Josephson effects. 30 ~ 32 ) 

Now, let us show a spatial modulation of the mean- field, i.e., the superconducting 
gap in the SNS junction with spatial size N x x N y = 64 X 32. The normal conducting 
part is inside the region {(x,y);27 < x < 37, 1 < y < N y }. We set Vij = inside 
this region. Otherwise, = —2.2t5ij. The chemical potential /i = — 1.5t and the 
temperature T = 0. The hopping between the superconducting and normal regions 
t' = 0.8t. We note that the averaged gap amplitude outside the normal region (i.e., in 
a uniform superconductor) is 0.198 in the present parameter set. It indicates that the 
coherence length is around 5 sites (Ref. 3). The mean-field superconducting-gap (q^q^) 
distribution is shown in Figs. 2(a) and (b). In this case, the self-consistent iteration 
number is 20, which is confirmed to be enough through its convergence check. The 
mean-field takes a non-zero value even in the normal region. We find that the spatial 
modulation of the mean-field is well characterized by the length scale compatible with 
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Fig. 3. (Color online) Temperature dependence of the gap amplitude in the two-dimensional s-wave 
superconductor. 



the coherence length (~ 5 sites). 



3.3 Temperature Dependence of Gap Amplitude 

We demonstrate temperature dependence of the superconducting gap amplitude. 
We consider a two-dimensional N x N square-plate s-wave superconductor. We set 
Vij = —2.2t5ij, the chemical potential fi = — 1.5t, and the system size N x N = 28 x 28. 
The 30 and 300 times iterations are adopted to evaluate the temperature dependent su- 
perconducting gap, respectively. As shown in Fig. 3, we successfully obtain the temper- 
ature dependence of the gap amplitude, 33 ) although the convergence tendency depends 
on the temperature range. These obtained temperature dependences are fitted by an an- 
alytical formula of the BCS gap function A(T) = A(0) tanh(1.74 v /A(0)/(1.764T) - 1). 
The slow convergence appears around T c , since the superconducting gap almost van- 
ishes near the critical point. However, this is not a fault of the expansion scheme. It is 
known that the diagonalization scheme also shows similar tendency. A calculation with 
300 times iterations takes about 3 minutes when using a parallel cluster with 112 cores. 

4. Technical remarks 

In this section, we refer to technical remarks to describe the advantages of the present 
polynomial expansion scheme. These points are useful when actually performing large- 
scale numerical calculations. 
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4 . 1 Matrix- Vector Product 

The calculations of p T q n m. Eq. (24c) is necessary to evaluate mean-fields or density 
of states. Since the target Hamiltonian is generally a sparse matrix, we can choose a fast 
algorithm optimized for the sparse-matrix vector product among several suggested ones. 
We confirm that Compressed Row Storage (CRS) format, which is one of the typical 
storing-schemes for sparse matrices, 34 ) is quite useful for the present calculations. The 
CRS format puts the subsequent nonzeros of the matrix row in contiguous memory 
locations. The algorithm of the matrix multiply operation y = Ax is easy for the 
program-coding. This algorithm is efficient on scalar processors since it has unit stride 
access. For details of the algorithm, see Ref. 34. 

4-2 Energy cut-off' 

In the present polynomial expansion scheme, the values of E max and E m \ n are initially 
demanded. A frequently-used way to do so is Lanczos 35 ) and its relative algorithm, as 
pointed out in the previous paper. 6 ) However, we mention that we do not need any 
exact calculations. Namely, what we need is just their approximate values. Tanaka et 
al. pointed out that calculation results are insensitive to the choice of these parameters. 
Therefore, we simply employ the band width as the upper and lower bounds of the 
energy range as E max = lOt — fi and E min = —lOt + fi, or E max = 20t — /i and E m i n = 
—20t + /i. In the previous section, we confirmed that this choice is effective. 

4-3 Remarks in Mean- field Calculations 

We discuss convergence tendency of the mean-fields. Equations (21a) and (21b) 
as a function of n are characterized by p T q n and T n , since w n is a constant for the 
Chebyshev polynomial for n > 1. With the use of the unitary matrix U diagonalizing 
the Hamiltonian ft, we obtain a recurrence formula with respect to x n = U^q n as 

[Xn+ll, = 2£ 7 [jc n ] 7 - [x n _i] 7 . (26) 

It indicates that 

2N 2N 

P T Qn = P^ U ^I COS [ n arCC0S (^)] • ( 27 ) 

H=l 7=1 

Therefore, p T q n is an oscillating function with respect to n. This is why the conversion 
of the local density of state (12) is slow, as Covaci et al. reported. 6 ) However, the 
conversion of the mean-fields is faster than that of the local density of states. As shown 
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Fig. 4. Schematic figure of N elements of the N x N Green's function. 



in Eq. (25b) at zero temperature, the n-dependence of T n is written as 

T n oc -. 

n 



(28) 



Then, p T q n T n shows a damped oscillating behavior on n. In the present calculations, 
we set n c = 1000, which is rather smaller than the value used in Ref. 6. On the one 
hand, we confirm that the value is enough for calculating the mean-fields. On the other 
hand, we find that the calculation of the local density of states requires larger n c than 
that of the mean-fields. When one tries to use the Chebyshev polynomial expansion 
scheme, one should choose a proper value of n c depending on the calculation target. 
However, the first choice is n c = 1000. Then, its tuning (e.g., n c = 500, 2000, 4000, 
etc.) should be done. One can easily try the various values since the computational 
time is proportional to the cutoff n c . Typically, the validity of the adopted value of n c 
is checked by investigating whether the resultant DOS has an oscillation behavior (i.e., 
Gibbs oscillation). 

4-4 Origin of Calculation Efficiency 

We mention a reason why the polynomial expansion scheme is much more efficient 
than full-diagonalization. Principally, the N xN matrix Green's function is constructed 
by all eigenvectors of "H. The full diagonalization directly calculates them. On the other 
hand, practically, all what we have to obtain to solve the BdG equations with the mean- 
fields is diagonal elements of the N x N Green's function as shown in Fig. 4. This 
indicates that the calculation cost can be considerably reduced compared to the full 
diagonalization. This is an origin of excellent efficiency of the polynomial expansion 
scheme. 

Here, let us evaluate computational costs in the present self-consistent calculations. 
We measure the elapsed time from making the Hamiltonian matrix with CSR format 
to finishing the 20 times self-consistent iterations in the M x M square lattice s-wave 
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Fig. 5. (Color online) System-size dependence of the self-consistent calculation in the M x M square 
lattice s-wave superconductor at zero temperature. The number of the iteration steps is 20. The system 
size denotes M x M. 



superconductor at zero temperature. For the measurement, we use a supercomputing 
system PRIMERGY BX900 in Japan Atomic Energy Agency. As shown in Fig. 5, the 
elapsed time of the self-consistent calculation grows in 0(N 2 ) manner with increasing 
the system size N(= M x M). The tendency is kept from 32 to 4096 cores. In fact, the 
computational cost is theoretically estimated to be 0(N 2 ) according to Eq. (21). Here, 
we mention that, although the sparse-matrix multiply operation is O(N), the mean field 
calculation on all sites requires an extra cost represented as O(N) x N. The previous 
report 6 -* claimed that the cost is O(N). However, it is practically 0(N 2 ) when includ- 
ing the mean field calculation. In contrast, the full diagonalization scheme inevitably 
demands 0(N 3 ) cost in the core part of the calculation. This is a big advantage of the 
polynomial expansion scheme. Furthermore, we focus on the speed of a self-consistent 
calculation by the present scheme. For example, a full calculation on 256 x 256 (2 16 ) 
lattice system, whose matrix dimension size is 131077, takes about 5 hours (~ 2 14 sec) 
for 20 iterations when using 1024 CPU cores. If we execute 20 times diagonalizations, 
then the elapsed time exceeds much over 5 hours on the same number of cores. The 
difference in CPU time becomes much more remarkable as the system size grows into 
nano to meso-scales. 

5. Conclusion 

In conclusion, we presented a full formulation of the polynomial expansion scheme 
for generic fermionic quadratic Hamiltonians with mean fields. The spectral density in 
Green's function was expanded by a set of polynomials, and the calculation scheme 
of the mean fields was explicitly described. The scheme was actually implemented to 
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solve the BdG equations for inhomogeneous superconductors. We demonstrated three 
examples of nano-scale self-consistent calculations for inhomegeneous superconductors, 
whose targets are a vortex lattice under randomly distributed impurities, proximity 
induced superconducting gap in an SNS junction, and temperature dependence of the 
gap amplitude in nano-square 2D plate superconductor. In all these calculations, we 
confirmed its high numerical efficiency We presented technical remarks to take ad- 
vantage of the polynomial expansion scheme. These practical remarks are useful when 
actually performing large-scale numerical calculations. The CPU cost scales with the 
system size N in the manner 0(N 2 ) in contrast to direct full diagonalization 0(N 3 ), 
and its relationship becomes much more crucial in larger systems. We claim that our 
scheme widely expands the calculation range and makes it possible to study meso-scale 
superconducting phenomena. 
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